package fem2.element.crack_sda2;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import math2.MathUtilities;
import fem2.Debugger;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;
import fem2.enu.EchoLevelType;
import fem2.enu.MeshPartType;
import fem2.enu.StressType;
import fem2.material.cohesive_law2.CohesiveLaw;

/**
 * CST element using Symmetric formulation and enforce crack line parallel with
 * element edge
 * 
 * @author hbui
 * 
 */
public class SOSConstantJumpCSTElement extends ConstantJumpSDAElement {

	public SOSConstantJumpCSTElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl) {
		super(mp, mm, cl);
		setCrackMeshPartType(MeshPartType.LINE2);
		setIntegrationOrderOnCrack(1);
	}

	/**
	 * SDA-SOS element using different stress type to evaluate crack
	 * 
	 * @param mp
	 * @param mm
	 * @param cl
	 * @param stressType
	 */
	public SOSConstantJumpCSTElement(MeshPart mp, MaterialModel mm, CohesiveLaw cl,
			StressType stressType) {
		this(mp, mm, cl);
		this.stressType = stressType;
	}

	@Override
	protected void computeEnhancedStrainOperator(MeshPart mp, IMatrix g, double... xi) {
		double c = -crackMeasure / elementMeasure;
		g.set(0, 0, c * normal[0]);
		g.set(1, 1, c * normal[1]);
		g.set(2, 0, c * normal[1]);
		g.set(2, 1, c * normal[0]);
		IMatrix J = computeLocalProjectionOnCrack();
		IMatrix tmp = (IMatrix) g.clone();
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, tmp, BLAM.NO_TRANSPOSE, J, 0.0, g);
	}

	@Override
	protected void computeProjectionOperator(MeshPart mp, IMatrix g, double... xi) {
		// computeEnhancedStrainOperator(mp, g, xi);
		double c = -crackMeasure / elementMeasure;
		g.set(0, 0, c * normal[0]);
		g.set(1, 1, c * normal[1]);
		g.set(2, 0, c * normal[1]);
		g.set(2, 1, c * normal[0]);
	}

	@Override
	protected boolean checkOnsetCriteria() {
		IVector sigma = getStress(0, stressType);

		double oxx = sigma.get(0);
		double oyy = sigma.get(1);
		double oxy = sigma.get(2);
		double t1 = (oxx + oyy) / 2;
		double t2 = Math.sqrt(Math.pow((oxx - oyy) / 2, 2) + Math.pow(oxy, 2));
		double o1 = t1 + t2;

		if (o1 < getCohesiveLaw().getTensileStrength()) {
			return false;
		} else {
			MeshPart mp = getMeshPart();
			int n = mp.getBasis().getDimension();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			alpha = new ArrayVector(dim);

			/*
			 * realize crack normal
			 */
			double n1, n2;
			if (Math.abs(oxy) < MathUtilities.tol14) {
				if (oxx > oyy) {
					n1 = 1;
					n2 = 0;
				} else {
					n1 = 0;
					n2 = 1;
				}
			} else {
				double g = Math.sqrt(Math.pow(oxx - o1, 2) + Math.pow(oxy, 2));
				n1 = -oxy / g;
				n2 = (oxx - o1) / g;
			}

			/*
			 * compute B operator
			 */
			Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			double[] xi = getIntegrationPoint(0);
			computeB(mp, B, xi);

			/*
			 * search the gradient which make smallest angle with crack normal
			 */
			double cos = 0.0;
			int posNode = 0;
			for (int i = 0; i < 3; i++) {
				double gradx = B.get(0, i * dim);
				double grady = B.get(1, 1 + i * dim);
				double normgrad = Math.sqrt(gradx * gradx + grady * grady);
				gradx /= normgrad;
				grady /= normgrad;
				double tmp = Math.abs(n1 * gradx + n2 * grady);
				if (tmp > cos) {
					cos = tmp;
					posNode = i;
					normal[0] = gradx;
					normal[1] = grady;
					/*
					 * compute crack length
					 */
					crackMeasure = elementMeasure * normgrad;
				}
			}
			tangent1[0] = -normal[1];
			tangent1[1] = normal[0];

			/*
			 * compute crack boundary
			 */
			int[] negNodes;
			if (posNode == 0) {
				negNodes = new int[] { 1, 2 };
			} else if (posNode == 1) {
				negNodes = new int[] { 0, 2 };
			} else if (posNode == 2) {
				negNodes = new int[] { 0, 1 };
			} else {
				throw new Error("positive node is wrong");
			}
			Node[] nodes = mp.getNodes();
			double[] X1 = new double[2];
			double[] X2 = new double[2];
			double[] posNodePos = nodes[posNode].getPosition();
			double[] negNode1Pos = nodes[negNodes[0]].getPosition();
			double[] negNode2Pos = nodes[negNodes[1]].getPosition();
			X1[0] = (posNodePos[0] + negNode1Pos[0]) / 2;
			X1[1] = (posNodePos[1] + negNode1Pos[1]) / 2;
			X2[0] = (posNodePos[0] + negNode2Pos[0]) / 2;
			X2[1] = (posNodePos[1] + negNode2Pos[1]) / 2;
			crackBoundary.add(X1);
			crackBoundary.add(X2);

			if (getEchoLevel().contain(EchoLevelType.AUX4)) {
				Debugger.watch("crack length = ", crackMeasure);
				Debugger.watch("crack normal = ", normal);
				Debugger.watch("element area = ", elementMeasure);
			}

			return true;
		}
	}

}
